%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%   原始数据  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 空气密度p
p = 1.25; 

% 爆炸时刻data_t以及对应时刻的蘑菇云半径data_r
data_t = 0.001 * [0.1 0.24 0.38 0.52 0.66 0.8 0.94 1.08 1.22 1.36 1.50 1.65 1.79 1.93 3.26 3.53 3.8 4.07 4.34 4.61 15 25 34 53 62];
data_r = [11.1 19.9 25.4 28.8 31.9 34.2 36.3 38.9 41 42.8 44.4 46 46.9 48.7 59 61.1 62.9 64.3 65.6 97.3 106.5 130 145 175 185];

% 蘑菇云半径与对应时刻的数据
figure(1);
plot(data_t, data_r, '*');
xlabel('t');
ylabel('r');
title('蘑菇云半径与对应时刻的数据');
saveas(1, '蘑菇云半径与对应时刻的数据.jpg');

% 计算ln(t)、ln(r)的值
data_t_log = log(data_t);
data_r_log = log(data_r);

% 使用一次多项式拟合ln(t)、ln(r)的关系
cof = polyfit(data_t_log, data_r_log, 1);
b = cof(1);
a = cof(2);

% 拟合数据并作图
figure(2);
plot(data_t_log, data_r_log, 'o', data_t_log, b * data_t_log + a, 'g-');
title('使用一次多项式拟合ln(r) = ln(a) + b*ln(t)的数据拟合图');
xlabel('ln t');
ylabel('ln r');
saveas(2, '使用一次多项式拟合的数据拟合图.jpg');

% 使用零次多项式拟合ln(E/P)、(5*ln(r)-2*ln(r))的关系
step = (1:1:25)';
data_t_r_log = polyfit(step, 5 * data_r_log - 2 * data_t_log, 0);
data_t_r_log_val = polyval(data_t_r_log, step);

% 拟合数据并作图
figure(3);
plot(step, 5 * data_r_log - 2 * data_t_log, 'o', step, data_t_r_log_val, 'r-');
title('使用零次多项式拟合(5*ln(r)-2*ln(r)) = ln(E/P)的数据拟合图');
xlabel('ln(E/P)');
ylabel('5*ln(r)-2*ln(r)');
saveas(3, '使用零次多项式拟合的数据拟合图.jpg');

% 根据E = p * e^c公式计算爆炸能量E
E = p * exp(data_t_r_log);
E_j = E / (4.184 * 10^(12));

fprintf("使用一次多项式拟合“ln(r) = ln(a) + b ln(t)”得到的系数为： %f\t%f\n", b, a)
fprintf("使用零次多项式拟合“(5*ln(r)-2*ln(r)) = ln(E/P)”得到的系数为： %f\n", data_t_r_log)
fprintf("通过公式E = p * e^c算出E的近似解为：%fJ\n", E)
fprintf("根据1kt = 4.184 * 10^(12) 转化公式，得到爆炸能量约等于： %fkt\n", E_j)
